//////////////////////////////////////////////////////////////////////////////////
//                                                                              //
//  This file is part of the buola project (https://code.google.com/p/buola/).  //
//                                                                              //
//  Copyright(c) 2007-2012 Xavi Gratal                                          //
//  gratal AT gmail DOT com                                                     //
//                                                                              //
//  Buola is free software: you can redistribute it and/or modify               //
//  it under the terms of the GNU General Public License as published by        //
//  the Free Software Foundation, either version 3 of the License, or           //
//  (at your option) any later version.                                         //
//                                                                              //
//  Buola is distributed in the hope that it will be useful,                    //
//  but WITHOUT ANY WARRANTY; without even the implied warranty of              //
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               //
//  GNU General Public License for more details.                                //
//                                                                              //
//  You should have received a copy of the GNU General Public License           //
//  along with buola.  If not, see <http://www.gnu.org/licenses/>.              //
//                                                                              //
//////////////////////////////////////////////////////////////////////////////////

#include <buola/image/rasterizer/cspanrasterizer.h>

namespace buola { namespace img {

static inline double gaussian(double r, double x, double y)
{
    return 1/sqrt(2*M_PI*r)*exp(-(x*x+y*y)/(2*r*r));
}

static img::CImage_gray64f sData;
static img::CImage_gray8 sCorners[26];
static img::CImage_gray8 sTops[26];

static void generate_data(double pR)
{
    int lSize=((int)ceil((pR*3))+1)&~1;
    int lCenter=lSize/2;

    sData.Recreate({lSize,lSize});

    double lG;

    double lT = 0.0;
    for(int y=0;y<lSize;y++)
    {
        for(int x=0;x<lSize;x++)
        {
            double lG=gaussian(pR,double(x-lCenter),double(y-lCenter));
            lT+=lG;
            sData(x,y)=lG;
        }
    }

    for(int y=0;y<lSize;y++)
        for(int x=0;x<lSize;x++)
            sData(x,y)/=lT;
}

/*
 * A picture will help
 *
 *      -center   0             width  width+center
 *  -center +-----+-------------------+-----+
 *          |     |                   |     |
 *          |     |                   |     |
 *      0   +-----+-------------------+-----+
 *          |     |                   |     |
 *          |     |                   |     |
 *          |     |                   |     |
 *   height +-----+-------------------+-----+
 *          |     |                   |     |
 * height+  |     |                   |     |
 *  center  +-----+-------------------+-----+
 */

static unsigned char sum(double pOpacity,int pX,int pY,int pW,int pH)
{
    int lFX,lFY;
    int lCenter=sData.Width()/2;
    int lFXStart,lFXEnd;
    int lFYStart,lFYEnd;
    double lV;

    /*
     * Compute set of filter values which are "in range",
     * that's the set with:
     *      0 <= x + (fx-center) && x + (fx-center) < width &&
     *  0 <= y + (fy-center) && y + (fy-center) < height
     *
     *  0 <= x + (fx - center)      x + fx - center < width
     *  center - x <= fx        fx < width + center - x
     */

    lFXStart=lCenter-pX;
    if (lFXStart < 0)
        lFXStart = 0;
    lFXEnd = pW + lCenter - pX;
    if (lFXEnd > sData.Width())
        lFXEnd = sData.Width();

    lFYStart = lCenter - pY;
    if (lFYStart < 0)
        lFYStart = 0;
    lFYEnd = pH + lCenter - pY;
    if (lFYEnd > sData.Width())
        lFYEnd = sData.Width();

    auto lIt=
    lGLine = lGLine + lFYStart * sData.Width() + lFXStart;

    lV = 0;
    for (lFY = lFYStart; lFY < lFYEnd; lFY++)
    {
        auto lIt=sData.XIterator(lFXStart,lFY);

        for (lFX=lFXStart;lFX<lFXEnd;lFX++)
            lV+=*lIt++;
    }
    if(lV>1)
        lV=1;

    return ((unsigned char) (lV * pOpacity * 255.0));
}

/* precompute shadow corners and sides to save time for large windows */
static void presum()
{
    int lCenter=sData.Width()/2;

    for(int i=0;i<26;i++)
    {
        sCorners[i].Recreate({sData.Width()+1,sData.Height()+1});
        sTops[i].Recreate({sData.Width()+1,1});
    }

    for (int x=0;x<=sData.Width();x++)
    {
        sTops[25](x,0)=sum(1.0,x-lCenter,lCenter,2*sData.Width(),2*sData.Width());
        for(int lOpacity=0;lOpacity<25;lOpacity++)
            sTops[lOpacity](x,0)=sTops[25](x,0)*lOpacity/25.0;
        for(int y=0;y<=x,y++)
        {
            sCorners[25](x,y)=sCorners[25](y,x)=sum(1.0,x-lCenter,y-lCenter,2*sData.Width(),2*sData.Width());
            for(int lOpacity=0;lOpacity<25;lOpacity++)
                sCorners[lOpacity](x,y)=sCorners[lOpacity](y,x)=sCorners[25](y,x)*lOpacity/25.0;
        }
    }
}

CImage_gray8 gaussian_map(double pRadius,double pOpacity,const CSize_i &pSize)
{
    if(sData.Null())
    {
        generate_data(10.0);
        presum();
    }
    
    /*
     * Build the gaussian in sections
     */
    int             lCenter=sData.Width()/2;
    int             lSWidth=pSize.x+sData.Width();
    int             lSHeight=pSize.y+sData.Height();
    int             lOpacityInt=(int)(pOpacity * 25);

    img::CImage_gray8 lImage({lSWidth,lSHeight},img::pixel((uint8_t)pOpacity*255.0));

    /*
     * corners
     */
    int lYLimit = sData.Width();
    if (lYLimit>lSHeight/2)
        lYLimit=(lSHeight+1)/2;
    int lXLimit = mSize;
    if (lXLimit>lSWidth/2)
        lXLimit=(lSWidth+1)/2;

    for(int y=0;y<lYLimit;y++)
    {
        for(int x=0;x<lXLimit;x++)
        {
            uint8_t lD;
            if(lXLimit==sData.Width()&&lYLimit==sData.Height())
                lD=sCorners[lOpacityInt](x,y);
            else
                lD=sum(pOpacity,x-lCenter,y-lCenter,pSize.x,pSize.y);
            lImage(x,y)=lD;
            lImage(x,lSHeight-1-y)=lD;
            lImage(lSWidth-1-x,lSHeight-1-y)=lD;
            lImage(lSWidth.x-1-x,y)=lD;
        }
    }
    /*
     * top/bottom
     */
    for(int y=0;y<lYLimit;y++)
    {
        uint8_t lD;
        if (lYLimit==sData.Width())
            lD=sTops[lOpacityInt](y,0);
        else
            lD=sum(pOpacity,lCenter,y-lCenter,pSize.x,pSize.y);
        for(int x=sData.Width();x<lSWidth-sData.Width();x++)
        {
            lImage(x,y)=lD;
            lImage(x,lSHeight-y-1)=lD;
        }
    }

    /*
     * sides
     */

    for(int x=0;x<lXLimit;x++)
    {
        uint8_t lD;
        if (lXLimit==sData.Width())
            lD=sTops[lOpacityInt](x,0);
        else
            lD=sum(pOpacity,x-lCenter,lCenter,pSize.x,pSize.y);
        for(int y=sData.Height();y<lSHeight-sData.Width();y++)
        {
            lImage(x,y)=lD;
            lImage(lSWidth-x-1,y)=lD;
        }
    }

    return lImage;
}

/*namespace img*/ } /*namespace buola*/ }
